Enhancing Ionic Conductivity of Bulk Single Crystal Yttria-Stabilized Zirconia by 

Tailoring Dopant Distribution 
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We present an ab-initio based kinetic Monte Carlo model for ionic conductivity in single crystal 
yttria-stabilized zirconia. Ionic interactions are taken into account by combining density functional 
theory calculations and the cluster expansion method and are found to be essential in reproducing 
the effective activation energy observed in experiments. The model predicts that the effective energy 
barrier can be reduced by 0.15-0.25eV by arranging the dopant ions into a super-lattice. 

PACS numbers: 82.20.Wt, 66.30.-h, 82.47.Ed, 82.45.Gj 



Yttria-stablized zirconia (YSZ) is a widely used elec- 
trolyte in solid oxide fuel cells (SOFC) and oxygen sen- 
sors because of its high ionic conductivity at high tem- 
peratures [l|. Driven by the need to reduce the operat- 
ing temperature of SOFC, much of the current research 
effort focuses on the design of new solid electrolyte ma- 
terials with significantly enhanced ionic conductivity at 
intermediate temperatures [1, 01 • The presence of free 
surfaces in nanoscale thin films and interfaces in hetero- 
epitaxial structures has been found to enhance ionic con- 
ductivity However, the effect of dopant distribu- 
tion on the ionic conductivity of bulk single-phase elec- 
trolytes has largely remained unexplored, despite the fun- 
damental importance of dopant-vacancy interaction in 
ionic transport [l^, [HI and the possibility of tailoring 
dopant distributions by novel deposition techniques 12 1 . 

Atomistic simulations have the promise to become 
a useful design tool for new electrolyte materials, by 
predicting the ionic conductivity of candidate struc- 
tures and elucidating the fundamental transport mech- 
anisms [III m " 12 1 • Unfortunately they are still limited 
in their length and time scales. For example, to accu- 
rately describe the long-range ionic interactions in YSZ 
requires density functional theory (DFT) models with 
relatively large supercells. The high computational cost 
limits the time scale of ab initio molecular dynamics sim- 
ulations to picoseconds [1^. Hence, a major challenge 
at present is to construct a kinetic Monte Carlo (kMC) 
model that not only can access the macroscopic time 
scale [H, 17|, but also retains the accuracy of DFT mod- 
els in describing the ionic interactions. In the pioneering 
kMC model for YSZ [l^ , ionic interactions are ignored in 
the metastable states, i.e., all possible states are sampled 
with uniform probability. While it successfully predicts 
a maximum in the conductivity as a function of doping 
concentration, the predicted temperature dependence is 
significantly weaker than experiments, signaling the im- 
portance of ionic interactions. The lack of ionic inter- 
action also makes this model unsuitable to predict the 
effect of dopant distribution on ionic conductivity. 

In this letter, wc develop a kMC model for oxygen va- 
cancy diffusion in YSZ that faithfully captures the ionic 



interactions. DFT calculations with supercell sizes sig- 
nificantly larger than previous studies [if 



18[ are per- 



formed to accommodate long-range interactions, and the 
data are used to construct a cluster expansion (CE) 
model. KMC simulations using this model predicts an 
effective activation energy that agrees better with exper- 
iments than the non-interacting model. The kMC simu- 
lations further predict that the maximum conductivity is 
achieved when the Yttrium dopant ions are distributed 
as [100] lines and form a 2D rectangular super-lattice in 
the two other directions. The effective energy barrier in 
this structure is lower than the random distribution by 
0.15-0.25eV. 

Ionic conduction in YSZ occurs through oxygen an- 
ion diffusion by the vacancy mechanism. The ionic con- 
ductivity is the averaged effect of many vacancy jumps 
and can be predicted from a kMC simulation over a 
sufficiently long time. The degrees of freedom in our 
kMC model are the positions of the oxygen vacancies, 
which hop on the simple cubic anion sublattice of YSZ. 
At each kMC step, the probability rates of all vacancy 
jumps to their nearest neighbor positions are calculated 
by j = '^o cxp ^— j^^^ , where ks is Boltzmann's con- 
stant, T is temperature, and Eb is the activation energy 
barrier for each jump, vo is a trial frequency and is set 
to lO^'^ Hz At each step, only one event is selected 
based on the probability rates of all possible events [3| ■ 
From a long kMC simulation, the diffusion coefficient of 
the vacancies center-of-mass is computed from the mean- 
square-displacement. The ionic conductivity is then com- 
puted from the Nernst-Einstein relation 27 1. 

A fundamental input to the kMC simulation is the en- 
ergy barriers for vacancy jumps, Eh, which depends on 
the ionic configurations around the jumping vacancy. In 
a previous model [3|, E^ is assumed to depend only on 
the chemical species of the two cations closest to the 
jumping vacancy, as shown in Fig. [TJa). Because the 
energy barrier of every jump equals that of the reverse 
jump, one can show that all metastable states in this 
model must have identical energy. Hence we will refer to 
it as the non-interacting model. Experimental and com- 
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FIG. 1: (a) In the non-interacting model Ei, depends only on 
whether cations A and B are Zr or Y. Et — 0.58, 1.29 or 
1.86 eV if within pair A-B there are 0, 1, or 2 Y cations [l^ . 
Filled circles represent cations, open circle represents oxygen 
anion, and open square represents oxygen vacancy, (b) In 
our interacting model, the energy difference AE between the 
initial and final state is assumed to modify the energy profile 
of vacancy jump by superimposing a linear function. 



putational data have suggested that interactions play an 



important role in ionic conduction jlOl . To account 
for interactions, we use the kinetically resolved activation 
(KRA) model in which Eb = f{E^,AE), where 

E^ is the "kinetically resolved" barrier when the two 
metastable states happen to have identical energy, and 
AE is the energy difference between the two states. We 
take the energy barriers in the non-interacting model [3| 
as our E^. The function / is often approximated by 
E, 



AE/2 [20|, 121 



Here we use a slightly bet- 
ter approximation for function / by assuming that the 
energy landscape between the two metastable states has 
a sinusoidal shape when AE = 0, and is modified by a 
linear term when AE ^ [27|. Hence our task of spec- 
ifying the energy barrier Eb is reduced to an accurate 
description of the energy difference AE. 

Because the ionic interactions in YSZ are long ranged, 
they can only be captured accurately by DFT calcula- 
tions in relatively large supercells. It is impossible to per- 
form DFT calculations for all ionic configurations sam- 
pled by the kMC simulation. Instead, we use the cluster 
expansion method (CEM) to limit the necessary number 
of DFT calculations. In CEM, every metastable state in 
the kMC simulation can be uniquely mapped to a spin 
configuration {si\ of an Ising model |18| . The energy as 



a function of ionic configuration can be expressed by a 
cluster expansion [23) 



E{{s.,}) ^Y.'^Mis.}) 



(1) 



where Va is called effective cluster interaction (ECI) for 



cluster a, and 



s, is the cluster function in- 



volving all spin variables belonging to cluster a. Eq. ([T]) 
is used to compute the energy of any metastable state 
encountered by the kMC simulation, prior to which the 
ECI coefficients V^s are fitted to a DFT data set. 

To construct the ab initio database for the fitting, we 
performed DFT calculations for 140 randomly chosen 
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FIG. 2: (color online) The 140 data points are divided into 
two sets: set I contains 100 data points (x) and is used to fit 
CEM, set II contains 40 data points (o) and is not included 
in the fit. Energies predicted from CEM are plotted against 
energies computed from DFT using VASP, in units of eV per 
cation. The inset shows the effective binding energy between 
Y and oxygen vacancy predicted by CEM. 



ionic configurations within a 3[100] x 3[010] x 3[001] YSZ 
supercell, using the Vienna Ab-initio Simulation Pack- 
age (VASP) [23[. This supercell contains 108 cation sites 
(Zr or Y) and 216 anion sites (O or Vg) and is signifi- 
cantly larger than previous studies [H, [ill, in order to 
accurately account for long range Coulombic and elastic 
interactions. K-point sampling is limited to the F-point 
considering the large size of supercell. The volume and 
the shape of the supercell are allowed to relax together 
with the ionic positions. A high energy cut-off of 520 eV 
is used to avoid Pulay stress. Each ionic configuration 
takes ^ 10** CPU-hours to be fully relaxed. 

The number of clusters needs to be significantly trun- 
cated for robust fitting, otherwise the cluster expansion 
model can be overly-adapted to the fitted data set (23 |. 
First, we only keep the clusters in which any two spins 
are separated by less than 1.5ao, where oq is the lattice 
parameter of YSZ. Second, we only keep clusters that in- 
volve up to 3 spins. Accounting for the translational and 
rotational symmetries, 192 independent clusters survive 
this truncation. For further truncation, a Monte Carlo 
algorithm is used to select Uc clusters out of 192 pos sible 
clusters by minimizing the cross validation score jl8| . To 
measure the predictive power of CEM, we separate the 
DFT data into two sets. Set I contains 100 data points 
and are used to fit the ECI's. Set II contains 40 data 
points are used to benchmark the CEM's predictions. 
The root mean square difference between the DFT ener- 
gies and the CEM's predictions per cation in Set I and 
Set II are defined as the error of fitting and the error of 
prediction, respectively. 

These two errors have different dependence on ric- For 
example, when ric = 97, the error of fitting is 0.0006 eV 
while the error of prediction is 0.012 cV. The large differ- 
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FIG. 3: (a) Conductivity of YSZ single crystal predicted 
by kMC simulations at T = 1800 K as a function of dop- 
ing concentration, (b) Predicted temperature dependence of 
YSZ conductivity at 8-mol% doping concentration. Y dopant 
cations are assumed to be randomly distributed. Experimen- 
tal data is reproduced from [2^ . 



ence between the two errors means that the fitted CEM 
has entirely lost its predictive power if ric is too large. 
Only when ?ic < 9, both errors are the same, and de- 
crease with increasing Uc- Hence, in this work, the op- 
timal choice is 71c = 9 [111 j where both the error of fit- 
ting and the error of prediction is 0.005 eV, as shown in 
Fig. [2] This error is small enough for our kMC simula- 
tions and is smaller than a previous study |18j . To our 
knowledge this is the first time the predictive power of 
an interaction model for YSZ has been demonstrated by 
monitoring the error of prediction. Fig. [H^b) shows the 
effective binding energy between an Y ion and an oxygen 
vacancy predicted by the CEM model. The preference of 
the oxygen vacancy to the second nearest neighbor site of 
Y is clearly seen, consistent with previous experimental 
Upl- 



and theoretical works 

The fitted CEM allows us to compute the energy dif- 
ference AE between the two states before and after a 
vacancy jump, which modifies the energy barrier Eb of 
the jump. Using this energy barrier model, we per- 
formed kMC simulations in 3[100] x 3[010] x 3[GG1] YSZ 
supercells in which the doping concentration varies from 
5% to 13%. The ionic conductivity at each doping con- 
centration is computed by averaging over 40 randomly 
generated Y distributions. As shown in Fig. ^a.), the 
ionic conductivity (at 1800 K) is maximum at 8-mol%, 



consistent with earlier experimental |25l . I2a | and theo- 



retical results [3 113 ■ -f'ig- ISIb) is the Arrhenius plot 
of ionic conductivity at 8mol% doping concentration. 
The predicted activation energy is 0.74eV at high T and 
0.85cV at low T, in much better agreement with experi- 
ments (0.85-1.0 eV) [2^ than the non-interacting model 
(0.59 eV) The remaining difference with experi- 

ments may be due to the error in E^ taken from [l6| . 
The neglect of activation entropy in vacancy jumps does 
not affect the conclusion because it does not change the 
slope of the Arrhenius plot. 

The conductivity results shown in Fig. |3] are the av- 



eraged value over 40 random Y configurations, the stan- 
dard deviation over which is about 19% of the average. 
This indicates that the ionic conductivity is sensitive to 
the spatial distribution of Y cations and poses the ques- 
tion: what is the optimal Y distribution that maximizes 
ionic conductivity? To answer this question, we have per- 
formed kMC simulations for a variety of Y distributions, 
in which the Y cations are segregated into either spherical 
clusters. (001) layers, or [100] rods. The simulation re- 
sults suggest two design principles that ultimately guide 
us to the optimal Y distribution. 

When all Y cations in the supercell are segregated 
into a spherical cluster, the ionic conductivity is actually 
lower than the random distribution (by 27% at 1800 K), 
contraryto the prediction based on the non-interacting 
model [15[. This is because the interaction between Y 
cations and oxygen vacancies, as shown in Fig. [2l at- 
tracts the vacancies next to the cluster. To diffuse over 
long distances and contribute to the ionic conductivity, 
vacancies must detach from the Y clusters. This requires 
overcoming a binding energy of ^ 0.12 eV, which reduces 
the ionic conductivity. 

Here the reduction of ionic conductivity is caused by 
the increase of spatial variation of the potential energy. 
Based on this result, we can formulate design principle 
I: the optimal Y distribution should minimize the energy 
variation of the metastable states as oxygen vacancies 
jump along the conduction direction. In this work, we 
focus on conduction along the [100] (i.e. x) direction. 
A promising candidate structure is to have Y cations 
segregated into planes, so that each (001) cation layer 
is either completely filled by Y, or completely filled by 
Zr. Due to translational invariance, the potential energy 
from cation-vacancy interaction remains constant after 
an oxygen vacancy jumps in the [100] direction 28[. The 
layered structure can be fabricated using thin film depo- 
sition techniques such as PLD [l2| . 

Unfortunately, the layered structure also has a lower 
conductivity than the random distribution (by 59% at 
1800 K). This is surprising because one might expect en- 
hanced conductivity, as the non-interacting model would 
predict, due to the existence of Y-free channels. The re- 
duction in conductivity is caused by the segregation of 
oxygen vacancies to the two anionic layers immediately 
next to the Y layer. Because there is always a first nearest 
neighbor (Inn) Y, vacancy diffusion in this layer experi- 
ences a high energy barrier (with E^ = 1.29 eV). Given 
that vacancies prefer to be the second nearest neighbors 
(2nn) of Y, as shown in Fig. [51 it is somewhat surpris- 
ing that the vacancies segregate to the nearest anionic 
layer. This problem is resolved by noticing that when 
the vacancy becomes the Inn of two Y cations, it be- 
comes the 2nn to four Y cations on the same cation layer. 
This result motivates our design principle II: the optimal 
Y distribution should not induce a high oxygen vacancy 
density in the first nearest neighbor sites of any Y cations. 
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FIG. 4: (color) (a) Vacancy distribution in structure B at 
T = 500 K - blue and gray circles indicate Y and Zr, respec- 
tively. Darker color in the orange columns corresponds to 
higher vacancy concentration, (b) Temperature dependence 
of the conductivity at different Y distributions. Random Y 
distribution corresponds to a doping concentration of 8mol%. 

It turns out that the design principles I and II can be si- 
multaneously satisfied when Y cations are segregated into 
[100] lines. When Y lines are well separated from each 
other, kMC simulations show that the oxygen vacancy 
density is peaked at 2nn sites surrounding the Y lines. 
To enhance conductivity, we would like to pack more Y 
lines per unit volume in order to enhance the overall va- 
cancy density. But when the distance between Y lines 
becomes too small, oxygen vacancies start to occupy Inn 
sites around Y, deteriorating the ionic conductivity. 

We examined a variety of 2D structures formed by [100] 
Y-segregated lines. The optimal structure at 1800 K is 
a square lattice with periodicity 1.5ao in y and z direc- 
tions (structure A). Interestingly, the optimal structure 
at 500 K is different; it is a rectangular lattice with peri- 
odicity l.Soo and 2oo in y and z directions, respectively 
(structure B). In both structures, the vacancy concen- 
tration is high in fast diffusing channels {E'^ = 0.58 eV) 
away from Inn sites of Y, see Fig. IHa). 

Fig.mjb) plots the temperature dependence of conduc- 
tivity for structures A and B. Their activation energy is 
around ~ 0.6 cV, significantly lower than the random 
distribution. Compared with the random Y distribution 
(at 8mol%), the ionic conductivity of structure A is en- 
hanced by a factor of 1.35 at 1800 K, 11 at 500 K and 
86 at 300 K. For structure B, the enhancement factor is 
22 at 500 K and 532 at 300 K. This result provides a 
theoretical upper limit of ionic conductivity that can be 
achieved by rearranging dopants in YSZ. 

In summary, wc have developed an ab-initio based 
kMC model of the vacancy diffusion in bulk YSZ that 
accurately accounts for the ionic interactions. The pre- 
dicted ionic conductivity shows much better agreement 
with experiments regarding the temperature dependence 
than the non-interacting model. The model predicts 
strong dependence of ionic conductivity on the spatial 
distribution of dopant cations. The maximum conduc- 



tivity is reached when Y cations are arranged into a rect- 
angular superlattice of [100] lines. Fabrication of this 
structure is challenging, but may be feasible with novel 
deposition techniques. The method presented here can be 
easily applied to other solid electrolytes, in which the op- 
timal dopant microstructure may be different from that 
in YSZ and may be easier to synthesize. 
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